Comprehensive analysis of cuproptosis-related prognostic gene signature and tumor immune microenvironment in HCC

Background: Copper is an indispensable mineral element involved in many physiological metabolic processes. Cuproptosis is associated with a variety of cancer such as hepatocellular carcinoma (HCC). The objective of this study was to examine the relationships between the expression of cuproptosis-related genes (CRGs) and tumor characteristics, including prognosis and microenvironment of HCC. Methods: The differentially expressed genes (DEGs) between high and low CRGs expression groups in HCC samples were identified, and further were analyzed for functional enrichment analysis. Then, CRGs signature of HCC was constructed and analyzed utilizing LASSO and univariate and multivariate Cox regression analysis. Prognostic values of CRGs signature were evaluated by Kaplan-Meier analysis, independent prognostic analysis and nomograph. The expression of prognostic CRGs was verified by Real-time quantitative PCR (RT-qPCR) in HCC cell lines. In addition, the relationships between prognostic CRGs expression and the immune infiltration, tumor microenvironment, antitumor drugs response and m6A modifications were further explored using a series of algorithms in HCC. Finally, ceRNA regulatory network based on prognostic CRGs was constructed. Results: The DEGs between high and low CRG expression groups in HCC were mainly enriched in focal adhesion and extracellular matrix organization. Besides, we constructed a prognostic model that consists of CDKN2A, DLAT, DLST, GLS, and PDHA1 CRGs for predicting the survival likelihood of HCC patients. And the elevated expression of these five prognostic CRGs was substantially in HCC cell lines and associated with poor prognosis. Moreover, immune score and m6A gene expression were higher in the high CRG expression group of HCC patients. Furthermore, prognostic CRGs have higher mutation rates in HCC, and are significantly correlated with immune cell infiltration, tumor mutational burden, microsatellite instability, and anti-tumor drug sensitivity. Then, eight lncRNA-miRNA-mRNA regulatory axes that affected the progression of HCC were predicted. Conclusion: This study demonstrated that the CRGs signature could effectively evaluate prognosis, tumor immune microenvironment, immunotherapy response and predict lncRNA-miRNA-mRNA regulatory axes in HCC. These findings extend our knowledge of cuproptosis in HCC and may inform novel therapeutic strategies for HCC.


Introduction
Hepatocellular carcinoma (HCC) accounts for approximately 90% of liver cancers and is a common cause of death in cancer patients (Villanueva, 2019). The World Health Organization estimates that more than 1 million patients will die of liver cancer by 2030, according to annual projections (Villanueva, 2019). Liver cancer has a 5-year survival rate of 18% and is the second most deadly cancer after pancreatic cancer (Jemal et al., 2017). Effective treatments for advanced liver cancer are currently unavailable, and most patients can only use a palliative therapy such as chemoembolization. Therefore, it is necessary for us the development of new methods to effectively treat liver cancer is of utmost importance (Liver and Cancer, 2012;Llovet et al., 2016).
Previous studies showed that a variety of therapeutic methods involving programmed cell death such as pyroptosis and ferroptosis play an important role in inhibiting the occurrence and development of tumors (Fang et al., 2020;Li and Li, 2020;Tan et al., 2020). Among them, ferroptosis plays an important role in inhibiting tumor growth, and the enrichment of intracellular iron ions enhances the therapeutic effect of anti-cancer drugs . Pyroptosis also has effects on tumor cell proliferation, invasion, and metastasis, thereby affecting cancer prognosis (Fang et al., 2020;Tan et al., 2020). Copper is an indispensable trace element involved in various biological processes. The content of copper in the human body under normal circumstances is maintained at a steady state, but the content is modified under pathological conditions. For example, the content of copper in tumor tissue and serum in various cancers is significantly increased compared with the content in normal tissues (Yaman et al., 2007;Baharvand et al., 2014;Dragutinović et al., 2014;Kaba et al., 2015). The dysregulation of copper storage induces oxidative stress and cytotoxicity, while the effective regulation of copper content can affect cancer progression (Davis et al., 2020;Golonka and Vijay-Kumar, 2021). Based on this principle, a variety of copper chelators and copper ionophores have been developed for the replacement of therapies against tumors (Brady et al., 2014;Safi et al., 2014;Jemal et al., 2017;Roth and Decaens, 2017).
Recently, Tsvetkov et al. (2022) discovered a new copperdependent mode of controlled cell death that differs from known cell death forms: copper can bind directly to fatty acylation components of the tricarboxylic acid cycle (TCA), leading to proteotoxic stress and ultimately cell death. This process is called cuproptosis, and cuproptosis-related genes (CRGs) were also identified. The Cyclin-Dependent Kinase Inhibitor 2A/Multiple Tumor Suppressor 1 (CDKN2A) gene, also known as P16 gene, encodes multiple tumor suppressor 1 and belongs to the INK4 family (Serrano, 1997). Dihydrolipoamide S-acetyltransferase (DLAT) and dihydrolipoamide S-succinyltransferase (DLST) are components of the pyruvate dehydrogenase (PDH) complex. The oligomerization of DLAT is due to the integration of copper and fatty acylated proteins in the TCA (Tsvetkov et al., 2022). Glutaminase (GLS) is an enzyme that converts glutamine to glutamate. Glutamine generates ATP by entering the TCA cycle through the mitochondrial oxidative phosphorylation or by recycling reducing equivalents through lactate secretion (glutaminolysis), thereby contributing to the production of energy and building blocks in cancer cells (Le et al., 2012;Fendt et al., 2013;Zhang et al., 2019). Pyruvate dehydrogenase E1 subunit alpha 1 (PDHA1) is a key component of PDH, which converts pyruvate to acetyl-CoA (Ozden et al., 2014). Deficiency of PDHA1 leads to mitochondrial dysfunction and promotes glycolysis (Zhong et al., 2015). PDHA1 has an important biological significance in a variety of human tumors (Song et al., 2019). Seven genes (FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, and PDHB) conferred resistance tocuproptosis, while three genes (MTF1, GLS, and CDKN2A) sensitized the cells to cuproptosis (Tsvetkov et al., 2022). Copper importers (SLC31A1) and copper exporters (ATP7A andATP7B) are important factors in keeping the intracellular copper concentration (Lutsenko, 2010). Mutations in ATP7A and ATP7B genes were found to cause Menke's disease and Wilson's disease (Nevitt et al., 2012). NLRP3 may influence tumor immunity mainly by mediating tumor-infiltrating lymphocytes and macrophages (Ju et al., 2021). NFE2L2 pathway is an important modulator of cell homeostasis, associated with enhanced tumor growth and aggressiveness (Hellyer et al., 2021). However, only few studies on CRGs in HCC are currently available. The study of the role of CRGs in HCC might be potentially beneficial for the development of new targets in the treatment of HCC.
This study comprehensively investigated the molecular alterations and clinical relevance of CRGs in HCC. A total of 371 HCC patients were selected from The Cancer Genome Atlas (TCGA) database, classified into 2 cuproptosis-related subtypes by consensus clustering based on 19 CRGs, the signaling pathways of differential gene enrichment in the two subtypes were analyzed, and the differences in immune cell infiltration, important immune checkpoints, and m6A methylation-related genes between the two subtypes were measured. The CRGs associated with the prognosis of HCC patients were also analyzed by logrank test and univariate regression analysis. Gene and protein expression of prognostic CRGs in HCC cells or tissues were also assessed, as well as their association with anti-tumor drug sensitivity, m6A methylationrelated genes, and tumor microenvironment (TME). Prognostic models were developed to predict the overall survival (OS) and progression-free survival (PFS) in patients with HCC. A competing endogenous RNA (ceRNA) regulatory network was constructed to screen the lnRNA-miRNA-mRNA network that affected the prognosis of HCC patients. Our analysis highlighted the importance of CRGs in HCC development and laid the foundation for the therapeutic application of cuproptosis regulators in HCC.

Data sources and preprocessing
CRGs and clinical information of patients with HCC were obtained from the TCGA database (https://portal.gdc.cancer. Gov/) (Tomczak et al., 2015). A total of 371 HCC tissues and 50 paracancerous tissues were collected and used in this study. The data used in this study were standardized data per million transcripts, and the data distribution was close to normal distribution. R software (v4.0.3) was used to extract gene expression data, and a data matrix was constructed for further analysis.

Identification of molecular subgroups
Nineteen CRGs were first extracted from the TCGA expression matrix. According to the consistent clustering of these 19 genes, the R software package Consensus-ClusterPlus (v1.54.0) was used for the analysis of consistency, and the maximum number of clusters was 6 (Wilkerson and Hayes, 2010). The TCGA-LIHCREAD cases were then divided into 2 clusters based on the expression profiles of CRGs. This process was repeated 500 times to ensure the stability and reproducibility of the classification.

Identification and functional enrichment analysis of DEGs
DEGs between C1 and C2 subgroups were obtained using Limma package (version 3.40.2) in R software (Ritchie et al., 2015). Adjusted p values were analyzed in the TCGA database to correct for false positive results. "Adjusted p < 0.05 and log2 (fold change) > 1 or log2 (fold change) < −1″ was defined as the criteria for screening the differential expression of mRNA. GeneMANIA (http://www.genemania.org), which is a software that elucidates the relationship between genes and datasets by building gene interaction networks (Warde-Farley et al., 2010), was used to visualize the gene network of CRGs through physical interaction, co-expression, prediction, co-localization, and genetic interaction, and to evaluate their functions.
The GO function and the enrichment of KEGG pathways were analyzed by the "clusterProfiler" package in R software (Yu et al., 2012). Potential biological pathways were also identified by GSEA (http://software.broadinstitute.org/gsea/index.jsp) (Powers et al., 2018). According to the data of TCGA, the DEGs were divided into upregulated and downregulated groups. Each analysis performed 10,000 gene combination permutations to identify signaling pathways with significant changes; the genes were considered as enriched for meaningful signaling pathways when p.adjust <0.05 and FDR (false discovery rate) < 0.25. Statistical analysis and graphing were performed using the R package clusterProfiler (3.18.0).

Expression of CRGs and survival analysis
The expression of CRGs in 371 HCC tissues and 50 paired adjacent tissues was analyzed using the TCGA database. In addition, univariate Cox regression analysis was used to investigate the effect of CRGs on the prognosis of HCC. Kaplan-Meier curves, p values, and hazard ratios (HR) with 95% confidence intervals (CIs) were obtained by logrank test and univariate Cox regression. Five CRGs (CDKN2A, DLAT, DLST, GLS, and PDHA1) with higher hazard ratios were screened from the Cox regression analysis plot. The relationship between the prognostic CRGs and the OS rate in HCC patients was also analyzed, and the AUC (Area Under Curve) under the ROC curve was calculated.

Cell lines and culture conditions
The following HCC cell lines and hepatocyte were used in this study: Hep3B, Huh7, HepG2, and L02 (Chinese Academy of Medical Sciences, Beijing, China). All cell lines were maintained in Dulbecco's modified Eagle's medium (DMEM; Gibco, Grand Island, NY, United States) supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, United States of), 100 U/ml penicillin, and 100 U/ml streptomycin (Invitrogen, Carlsbad, CA, United States) and incubated at 37°C in a humidified atmosphere with 5% CO 2 .

Quantitative RT-PCR
First, TRIzol reagent (Invitrogen, CA, United States) was used to extract total RNA from the samples. Then, the total RNA concentrations were measured using a NanoDrop 2000c (Thermo Fisher Scientific, Inc.). Next, cDNA was synthesized by reverse

Development of the CRG prognostic model
A prognostic model was constructed using LASSO-Cox regression analysis based on the above CRGs associated with the prognosis of HCC patients. According to the results of the multivariate Cox regression analysis, the prognostic CRG risk score was calculated as follows: Risk score = i Cofficient(mRNAi) × Expression (mRNAi). Then, TCGA-LIHC patients were divided into low-risk and high-risk subgroups according to the average risk score, Kaplan-Meier analysis was used to compare the OS rates of the two subgroups, and timeROC analysis was performed to predict the accuracy of the model. Each variable (including the p-value, and HR with 95% CI) was displayed by univariate and multivariate cox regression analysis and using forest plots by the "forestplot" package. The "rms" package was used to build a Nomogram model for predicting 1, 3, and 5-year OS and PFS based on the results of multivariate cox proportional hazards analysis.

Tumor staging analysis of HCC
GEPIA is a newly developed interactive website for analyzing RNA sequencing expression data from 9736 tumors and 8587 normal samples from TCGA and Genotype-Tissue Expression databases, such as differential expression analysis of tumor/normal tissues, tumor type or pathological stage, patients' survival analysis, similar gene detection analysis, and dimensionality reduction analysis (Tang et al., 2017). The expression data of prognostic CRGs in different stages of HCC were obtained using the "Stage Plot" module in the GEPIA2 database (http://gepia2. cancer-pku.cn/#index). The UALCAN database (http://ualcan.path. uab.edu/) is an interactive portal containing TCGA RNA transcriptome data and clinical data from 31 cancer types. TISIDB (http://cis.hku.hk/TISIDB) is a portal for assessing tumors and the immune system that integrates multiple heterogeneous data types including genomics, transcriptomics, and clinical data from TCGA of 30 cancer types (Ru et al., 2019). The UALCAN and TISIDB databases were used to confirm the association between prognostic CRGs and the clinical stage of HCC patients. p < 0.05 was considered statistically significant.

Immunohistochemistry of prognostic CRGs in HCC
Immunohistochemical images of CRGs were collected from the Human Protein Atlas (https://www.proteinatlas.org) database to assess the differences in CRG expression between HCC and adjacent tissues (Sun et al., 2021).

Mutation analysis of CRGs
cBioPortal (http://www.cbioportal.org/) provides visualization tools for analyzing cancer genetic data. The genomic maps of CRGs were analyzed using cBioPortal based on the TCGA database to understand the mutation frequency in HCC (Gao et al., 2013). In this study, 372 LIHC samples was selected to explore the impact of prognostic CRGs on survival of HCC patients. In addition, the effects of CRG mutations on cancer signaling pathway expression, TMB, MSIsensor scores, and hypoxia-related scores (Ragnum Hypoxia Score, Buffa Hypoxia Score, and Winter Hypoxia Score) were examined. p < 0.05 was considered statistically significant.

Effects of prognostic CRGs on immune cell infiltration and immune checkpoint expression
The effects of prognostic CRGs on the abundance of infiltrating immune cells in tumors were analyzed by the TIMER (https:// cistrome.shinyapps.io/timer/) database . The detected immune cells included tumor purity, B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells. In addition, the infiltration of immune cell types was quantified by the R package "GSVA" single sample gene set enrichment analysis (ssGSEA) (Hanzelmann et al., 2013). Next, the Spearman method was used to explore the correlation of prognostic CRGs with the level of immune cell infiltration. Finally, the relationship between prognostic CRGs and checkpoints (CD274, CTLA4, and PDCD1) was examined.

TMB, microsatellite instability and drug sensitivity
The Spearman's method was used to analyze the correlation between prognostic CRGs and MSI in HCC. The chemotherapeutic response was predicted for each sample using the GDSC database (https://www.cancerrxgene.org/) (Yang et al., 2013). The R package pRRophetic was used to predict the half maximum inhibitory concentration (50% inhibition of the concentration, IC50) of the chemotherapeutic drugs, in which the IC50 of the sample was estimated by ridge regression. All parameters were set by default, using the batch effect of combat and the organization types of all, and summarizing the replicate gene expression as an average. Drug sensitivity and gene expression profiling data from cancer cell lines in the GDSC database were integrated in this study.
Frontiers in Genetics frontiersin.org

Single cell analysis
The effect of prognostic CRGs on the expression of single cell subsets in the TME was investigated using TISCH (http://tisch. comp-genomics.org/) (Sun et al., 2021). TISCH is a scRNA-seq database focused on the TME, providing detailed cell-type annotation at the single-cell level, which is beneficial for exploring the TME in different cancer types. In this dataset, three main cell types are present, such as immune cells, stromal cells, and malignant cells. In this study, the t-distributed stochastic neighborhood embedding (t-SNE) map of LIHC_GSE125449_

FIGURE 1
The flowchart of the present study.
Frontiers in Genetics frontiersin.org aPDL1aCTLA4 and the heatmap of LIHC_GSE125449_ aPDL1aCTLA4 were displayed through the TISCH database to demonstrate the effect of CRGs on the TME of HCC.

Identification and analysis of cuproptosisrelated gene clusters in HCC
The flowchart of the study is illustrated in Figure 1. A total of 371 HCC carcinoma samples were clustered in the TCGA database using consensus clustering T = to identify potential CRG clusters. All  Frontiers in Genetics frontiersin.org 07 tumor samples were divided into k (k = 2-6) distinct clusters based on the expression of 19 CRGs in HCC (Figures 2A,B). Then, the number of clusters was selected as 2 according to the cluster analysis results, indicating that the patients with HCC were accurately divided into two clusters (Cluster 1 had 252 patients and Cluster 2 had 119 patients) ( Figure 2C). Figure 2D shows the heat map of cuproptosis-related gene expression in different subgroups.

Differentially expressed genes and functional enrichment analysis
The study included 421 liver HCC samples consisting of 371 tumor samples and 50 adjacent normal samples. The DEGs identified between the C1 and C2 subtypes contained 5044 upregulated genes and 52 downregulated genes. Then, a volcano map ( Figure 3A) and a heat map ( Figure 3B) were constructed based on these DEGs. The results of GeneMANIA revealed that the functions of the significant genes co-expressed in this network (DLD, GLS, NFE2L2, DLST, LIAS, ATP7B, PDHA1, SLC31A1, FDX1, DBT, GCSH, LIPT1, CDKN2A, DLAT, NLRP3, ATP7A, PDHB, LIPT2, and MTF1) were related to several processes, including oxidoreductase complex, tricarboxylic acid cycle enzyme complex, cellular amino acid catabolic process, dihydrolipoyl dehydrogenase complex, acyl-CoA metabolic process, acetyl-CoA biosynthetic process and thioester metabolic process ( Figure 3C). The identified up and downregulated DEGs were then further subjected to GO and KEGG enrichment analysis. The results of biological process (BP) analysis showed that the DEGs were mainly enriched in extracellular matrix tissues, regulation of DNA metabolism, response to transforming growth factor β, activation of protein kinase activity, and cell cycle G1/S phase transition. The results of cellular component (CC) analysis showed that the DEGs were mainly enriched in focal adhesions, cell-matrix adhesion junctions, extracellular matrix components, collagen-containing extracellular matrix, and fibrous collagen trimers. The results of molecular function (MF) analysis showed that the DEGs were mainly enriched in several functions including cell adhesion molecule binding, extracellular matrix structural components, growth factor binding, protein serine/threonine kinase activity, and collagen binding. The KEGG analysis showed that the DEGs were mainly enriched in focal adhesion, extracellular matrixreceptor interaction, PI3K-Akt signaling pathway, TGF-β signaling pathway, HCC, cell cycle, NOD-like receptor signaling pathway, PD-L1 expression and the PD-1 checkpoint pathway in cancer ( Figure 3D). The results of Gene Set Enrichment Analysis (GESA) showed that CRGs were closely related to cancer pathways in HCC, including focal adhesions, cell cycle, T cell receptor signaling pathway, JAK STAT signaling pathway and MAPK signaling pathway. The activation of these pathways increased the risk of tumorigenesis and tumor progression ( Figure 3E).

Analysis of the correlation with immune infiltration and immune checkpoints
The TCGA-LIHC cohort was downloaded to explore the difference in immune cell infiltration between the C1 and C2 subtypes to explore the role of CRGs in the immune response in HCC. Six state-of-the-art algorithms were then integrated for a reliable immune score assessment. quantTIseq, EPIC, MCP-counter, CIBERSORT, xCell and TIMER showed the expression of CRGs was correlated with uncharacterized cell, Neutrophil, NK cell, Macrophage M1, T cell regulatory (Tregs), Macrophage M2, B cell ( Figures 4A,B) Figure S1E) respectively. Finally, the difference between the two subgroups in the expression of eight immune checkpoint-related genes was explored, and the results showed that CD274 (p < 0.01), HAVCR (p < 0.01), PDCD1LG2 (p < 0.01), TIGIT (p < 0.05) and SIGLEC15 (p < 0.05) were highly expressed in the C1 subgroup than in the C2 subgroup ( Figure 4C).

DEGs and prognostic models
The expression of 19 CRGs in HCC and normal tissues was investigated using the TCGA dataset; among them, 15 CRGs were up-or downregulated. The expression of ATP7A, SLC31A1, LIAS, LIPT1, LIPT2, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A, and DLST was upregulated, while NLRP3 and DBT were downregulated in cancer tissues than in the normal tissues ( Figure 5A). Most of the 19 CRGs in the HCC samples were positively correlated, while 6 pairs were negatively correlated, showing a strong correlation between them ( Figure 5B). Univariate Cox regression analysis was then performed to screen CRGs with a prognostic value ( Figure 5C). The mRNA expression of five prognostic CRGs was significantly upregulated (p < 0.001) in HCC cell lines (Hep3B, Huh7 and HepG2) compared to their expression in the hepatocyte cell line (L02) ( Figure 5D). Finally, five genes with a prognostic value were identified from the results. The Kaplan-Meier survival curves are shown in Figure 5E. HCC patients had a lower survival rate when CDKN2A (HR = 1.75, p = 0.002), DLAT (HR = 1.7208, p = 0.0024), DLST (HR = 1.5484, p = 0.0141), GLS (HR = 1.5001, p = 0.0222) and PDHA1 (HR = 1.4981, p = 0.0226) were highly expressed.

Construction of a prognostic CRG model
A prognostic gene model based on prognostic CRGs was constructed by LASSO Cox regression analysis (Figures 6A,B). The risk score was calculated as follows: (0.1242)*CDKN2A+ (0.2058)*DLAT + (0.0795)*GLS. HCC patients were divided into two groups according to the risk score distribution, survival status, and expression of GLS, DLAT, and CDKN2A ( Figures 6C,D). Patients had an increased risk of death and decreased survival Frontiers in Genetics frontiersin.org time as the risk scores increased ( Figure 6C). The Kaplan-Meier curves showed that HCC patients with high-risk scores had a lower probability of OS than patients with low-risk scores (median time = 3.4 years vs. 5.8 years, p = 0.00106) ( Figure 6D). The area under the ROC curve (AUC) in the 1-year, 3-year and 5-year ROC curves was 0.737, 0.646 and 0.634, respectively ( Figure 6E). The results of our  Frontiers in Genetics frontiersin.org 10 constructed cuproptosis-related risk profile showed a significant association with survival in patients with HCC.

Building of a predictive nomogram
Considering the clinicopathologic features and prognostic CRGs, a predictive nomogram was also built to predict the survival probability. The results of the univariate and multivariate analysis showed that CDKN2A expression and pT stage were independent factors affecting the prognosis of HCC patients ( Figures 7A,B; Supplementary Figure S2). The predictive nomogram suggested that the 3-year and 5-year OS rates and PFS rates were predicted enough accurately compared with an ideal model in the entire cohort (Figures 7C,D; Supplementary Figure S2).

Correlation between prognostic CRGs and pathologic stage in HCC
Based on the above analysis, the correlation between the expression of CRGs and the pathological stage of HCC was analyzed by Gene Expression Profiling Interactive Analysis (GEPIA) 2 database. The results showed that the expression of CDKN2A, DLAT, GLS and PDHA1 were significantly correlated with the pathological stage of HCC (p < 0.05), while the correlation between DLST and the pathological stage was not significant (Supplementary Figure S3). The analysis of UALCAN database revealed that the expression of prognostic CRGs in HCC patients in stages 1-3 was significantly upregulated compared with their expression in the normal group. The analysis of TISIDB database revealed that the expression of CDKN2A, DLST, and GLS in different stages and scores of HCC patients was significantly different, while DLAT and PDHA1 was slightly changed although not significantly (Supplementary Figure S3).

Pathological expression of CRGs in HCC tissues and normal livers
The protein expression of CRGs in 371 HCC and adjacent tissues was measured using the Human Protein Atlas to determine the difference in protein expression of the prognostic CRGs in HCC tissues and normal liver. The results of the immunohistochemical staining showed that the prognostic CRGs were moderately or highly expressed in HCC tissues but were low in paracancerous tissues ( Figures 8A,B).

TMB, microsatellite instability and drug sensitivity analysis
Microsatellite instability (MSI) and TMB can be used as indicators to predict the response of various tumors to immunotherapy (Rizzo et al., 2021). The above results showed that the expression of CRGs was closely related to tumor immune cell infiltration in HCC. Next, the association of prognostic CRGs with TMB and MSI in HCC was analyzed to determine whether CRGs can be considered as biomarkers for immunotherapy. The results showed that PDHA1 was positively Frontiers in Genetics frontiersin.org correlated with MSI (p = 0.044) ( Figure 11A). Furthermore, CDKN2A (p = 1.64 × 10 −04 ) was positively correlated with TMB, and GLS (p = 7.74 × 10 −05 ) and negatively correlated with TMB ( Figure 11B). Finally, the data on the gene expression profiles of cancer cell lines were integrated in the Genomics of Drug Sensitivity in Cancer (GDSC) database for drug sensitivity to fully explore the potential value of CDKN2A, DLAT, DLST, GLS, and PDHA1 genes as novel therapeutic targets. A Pearson correlation analysis of the data was performed, which showed that the expression of CDKN2A, DLAT, DLST, GLS, and PDHA1 was positively correlated with Saracatinib, Selumetinib, Serdemetan, Rucaparib, Roscovitine, Refametinib, Proteasome, Palbociclib, Motesanib, Lapatinib, Imatinib and Erlotinib, but negatively correlated with Tretinoin, Tipifarnib, Nilotinib, Navitoclax, Doramapimod, CCT018159 and AICA_ribonucleotide ( Figure 11C).

Correlation between CRGs and m6A methylation-related genes in HCC
It is well known that m6A methylation modification affects gene expression by regulating RNA metabolism and plays an important were found between the C1 and C2 subgroups. The results showed that except IGF2BP1, other m6A-related genes were significantly different between the two groups (p < 0.01, Figure 12A). The C1 group with higher expression of CRGs has also higher expression of m6A methylation-related genes. In addition, the correlation between prognostic CRGs and m6A-related genes was analyzed by the TCGA dataset, and the results showed a significant positive correlation between prognostic CRGs and m6A-related genes ( Figure 12B). The above results suggested that CRGs were closely related to m6A modification in HCC.

Single-cell RNA data analysis
The TME is composed of extracellular matrix, cancer associated fibroblasts (CAFs), myofibroblasts, immune cells and other factors. The prognostic CRGs were assessed at a single-cell level expression by Tumor Immunity Single Cell Center (TISCH, http://tisch.compgenomics.org/) to examine the relationship between prognostic CRGs and the infiltration of immune cells, stromal cells, and malignant cells. The HCC single-cell GSE dataset was explored (LIHC_GSE125449_aPDL1aCTLA4). T cells, B cells, plasma cells, Frontiers in Genetics frontiersin.org macrophages, endothelial cells, fibroblasts, malignant cells, and hepatic progenitors were annotated by single-cell RNA sequencing analysis ( Figure 13A). The results showed that CDKN2A, DLAT, DLST, GLS, and PDHA1 were correlated in single cell subsets such as immune cells, stromal cells, and malignant cells, and were significantly expressed in macrophages, fibroblasts, endothelial cells, and hepatic progenitor cells ( Figures 13B,C). CSFs and tumor-associated macrophages (TAMs) play an important role in the occurrence and development of cancer (Rizzo et al., 2021). Therefore, the association between prognostic CRGs and biomarkers associated with CAFs and TAMs was further explored. The results showed that prognostic CRGs were significantly and positively correlated with CAF markers such as FAP, PDPN, S100A8, S100A9, TGFB1, TGFB2, ACTA2, PALLD, TNC, and COL11A1 ( Figure 13D). Prognostic CRGs were also highly correlated with TAM-related markers such as CCL2, CD68, and IL 10 ( Figure 13E).

Prediction and validation of upstream key miRNAs
Next, the upstream regulatory miRNAs of prognostic CRGs were predicted using a comprehensive miRNA-related database.
Firstly, 27 pairs of CDKN2A-miRNAs, 63 pairs of DLAT-miRNAs and 172 pairs of GLS-miRNAs were obtained by the intersection of ENCORI and RNAInter databases. A total of 110 pairs of DLST-miRNAs and 39 pairs of PDHA1-miRNAs were obtained by the intersection of ENCORI databases and RNA22 databases ( Figure 14A). Then, according to the classical mechanism of miRNAs negatively regulating mRNA expression, a negative correlation between the mRNA and the predicted miRNA was expected. The ENCORI database was used to screen the correlation, prognosis and expression of these candidate miRNAs in HCC. Among these miRNA-mRNA interactions, the results showed that 1 pair of miRNA-CDKN2A, 6 pairs of miRNAs-DLAT, 3 pairs of miRNAs-DLST, 18 pairs of miRNAs-GLS, and 2 pairs of miRNAs-PDHA1 were significantly and negatively correlated ( Figure 14B; Supplementary Figure S6). Theoretically, miRNAs that bind to prognostic CRGs should be downregulated and indicate poor prognosis in HCC. Therefore, the prognostic role and expression of these potential miRNAs in HCC were further confirmed using the ENCORI database. The results showed that only the low expressed hsa-miR-125b-5p, hsa-miR-101-3p and hsa-miR-23c corresponded to a poor prognosis, and hsa-miR-125b-5p, hsa- Frontiers in Genetics frontiersin.org miR-101-3p and hsa-miR-23c were significantly lower in HCC tissues than in normal tissues ( Figure 14C). All these results indicated that CDKN2A-hsa-miR-125b-5p, GLS-hsa-miR-125b-5p, GLS-hsa-miR-101-3, GLS-hsa-miR-23c and PDHA1-hsa-miR-125b-5p might represent key pathways that mediated the occurrence and development of HCC and were related to the prognosis of patients.

FIGURE 12
Correlation of CRG expression and m6A-related genes in HCC. (A) Differential expression of m6A-related genes between the C1 and C2 subgroups in HCC. (B) Analysis of the relationship between prognostic CRGs and m6A-related genes by TCGA-LIHC cohort.
Frontiers in Genetics frontiersin.org 17 Prediction and validation of key lncRNAs binding to potential miRNAs Previous studies showed that lncRNAs bind to miRNAs and mediate the expression of target genes, thereby exerting biological roles (Militello et al., 2017;Zhang and Lou, 2020). Based on the above results, lncRNAs upstream of miRNAs were also predicted to construct the miRNA-lncRNA axis. The lncRNAs potentially binding to hsa-miR-101-3p, hsa-miR-125b-5p and hsa-miR-23c were predicted by the intersection of ENCORI and miRNet databases, and the results revealed 27 lncRNAs targeting hsa-miR-101-3p, 47 lncRNAs targeting hsa-miR-125b-5p and 56 lncRNAs targeting hsa-miR-23c ( Figure 15A). A miRNA-lncRNA regulatory network was established by Cytoscape software for a better visualization ( Figure 15B). According to the ceRNA hypothesis, lncRNAs can increase mRNA expression by competitively binding to miRNAs. Therefore, lncRNAs were negatively correlated with miRNAs or positively correlated with mRNAs. The correlation between lncRNAs and hsa-miR-101-3p, hsa-miR-125b-5p and hsa-miR-23c was detected by the ENCORI database, and the results showed that 14 lncRNAs were significantly associated with hsa-miR-101-3p and GLS, 9 lncRNAs were significantly associated with hsa-miR-125b-5p and CDKN2A, 10 lncRNAs were significantly associated with hsa-miR-125b-5p and GLS, 4 lncRNAs were significantly associated with hsa-miR-125b-5p and PDHA1, and 12 lncRNAs were significantly associated

Discussion
Cuproptosis is a novel cell death mechanism that is characterized by cell death induced by copper, which targets proteins of the fatty acylated tricarboxylic acid cycle (Tsvetkov et al., 2022). Wilson's disease is a condition caused by copper accumulation. Previous studies showed that humans or animal models with Wilson's disease have an increased incidence of HCC, indicating that the accumulation of copper promotes the malignant transformation, but the mechanism is not yet clear (Czlonkowska et al., 2018). CRGs play a crucial role in the development of kidney cancer (Bian et al., 2022). However, the role of CRGs in HCC has not yet been elucidated. Therefore, a bioinformatic analysis of public sequencing data of CRGs in HCC was performed in this work to gain a deeper understanding of CRGs expression, prognosis and their potential biological functions in HCC and RT-qPCR was used for experimental validation to guide future research.
The results showed that these 19 CRGs were mainly involved in extracellular matrix organization, regulation of DNA metabolism, cell-matrix adhesion junction, extracellular matrix components, cell adhesion molecule binding, collagen binding and other biological functions in HCC, and they were closely related to focal adhesion, extracellular matrix-receptor interaction, HCC, cell cycle, PD-L1 expression, PD-1 checkpoint, PI3K-Akt, TGF-β, NOD-like receptor and other signaling pathways. The results of GSEA enrichment analysis showed that the potential biological processes and pathways involved in CRGs in HCC included focal adhesion, cell cycle, T cell receptor, JAK STAT, MAPK and other signaling pathways. These pathways are closely related to tumorigenesis and tumor progression. The adhesion of cancer cells to the extracellular matrix promotes the resistance of cancer cells to chemotherapeutic drugs, thereby promoting the occurrence and development of tumors (Eke and Cordes, 2015). Signaling proteins at focal adhesions include kinases such as focal adhesion kinase, integrin-linked kinase, and phosphatases. Focal adhesion kinase is a cytoplasmic tyrosine kinase identified as a key mediator of integrin signaling (Seong et al., 2013). It plays an important role in tumor progression and metastasis by regulating cancer cell functions such as migration, invasion, and epithelial-mesenchymal transition, as well as affecting the pericancer microenvironment such as angiogenesis (Zhao and Guan, 2009). The focal adhesion kinase inhibitor BI853520 inhibits the proliferation, migration, invasion and epithelial-mesenchymal transition of cancer cells by affecting the PI3K/AKT/mTOR signaling pathway Li et al., 2021). Our results suggested that CRGs may affect the occurrence and progression of HCC by affecting focal adhesion and various cancer-related signaling pathways.
According to our results the expression of most CRGs in HCC tissues was significantly higher than that in the adjacent tissues. Quantitative RT-PCR results also confirmed that The expression of CDKN2A, DLAT, DLST, GLS, and PDHA1 was higher in various HCC cell lines than in the hepatocyte cell line (L02). The results of prognostic analysis showed that HCC patients with higher expression of CDKN2A, DLAT, DLST, GLS, and PDHA1 had shorter OS, indicating that prognostic CRGs promoted the progression of HCC. CDKN2A is highly expressed in various cancer tissues such as liver cancer and kidney cancer and plays an important prognostic role in many cancers (Ai et al., 2003;Zeng et al., 2018;Christodoulou et al., 2020;Ji et al., 2020;Xande et al., 2020). DLST is an E2 component of the α-ketoglutarate (αKG) dehydrogenase complex, which governs the entry of glutamine into the TCA for oxidative decarboxylation, thus promotes neuroblastoma aggression (Anderson et al., 2021). Over expression of glycolytic gene DLAT, which promoted glycolysis but suppressed acetyl-CoA production and enhanced the malignancy of non-small cell lung cancer (NSCLC) cells. Clinically, high expression of DLAT was positively associated with tumor size, poorer prognosis, and SUVmax values of 18F-FDG-PET/CT scans in patients with NSCLC (Chen et al., 2022). GLS the first enzyme in glutaminolysis, is overexpressed in ibrutinib-resistant mantle cell lymphoma (MCL) cells, and that its expression correlates well with elevated glutamine dependency and glutaminolysis. Targeting glutaminase is therapeutically effective in ibrutinib-resistant MCL (Li et al., 2022). PDHA1, a subunit of the pyruvate dehydrogenase complex (PDC), inhibits prostate cancer development in mouse and human xenograft tumor models by affecting lipid biosynthesis (Chen et al., 2018). Therefore, upregulated CRGs has an important biological significance in a variety of human tumors.
Univariate Cox regression and LASSO Cox regression analysis were subsequently used to analyze the five signature genes. The prognostic value of these five genes was then tested, and the results showed that HCC patients with high expression of these five CRGs had decreased survival rates. Moreover, the higher expression of GLS, DLAT, and CDKN2A increased the risk score, and the OS rate of patients in the high-risk group was significantly shorter than that of the patients in the low-risk group. In addition, CDKN2A was found as an independent risk factor, and a valid nomogram was constructed to predict the 1-, 3-, and 5-year survival rates of HCC patients, suggesting that CDKN2A played an important role in the occurrence and prognosis of HCC.
Frontiers in Genetics frontiersin.org Most malignancies are caused by somatic mutations within the cancer genome, and mutational signatures correlate with mRNA and protein expression in multiple cancer types (Weir et al., 2004). Functional alterations due to somatic mutations in cancer genomes are critical for identifying driver mutations and developing molecularly targeted therapies (Jia and Zhao, 2017). In this study, the analysis of the cBioPortal database showed that prognostic CRGs had a high frequency of gene alterations in HCC patients, and the survival rate of HCC patients with mutations in prognostic CRGs was significantly reduced compared with those without mutations. In other words, the mutation of prognostic CRGs in HCC accelerated the progression of HCC, providing some directions for the development of targeted drugs in the treatment of HCC in the future.
The interaction of the tumor and the immune system also plays an important role in the occurrence, development, and treatment of cancer (Kong et al., 2020;Xiao et al., 2020). Immune cells are an important part of the TME, and innate immune cells (including macrophages, neutrophils, dendritic cells, and natural killer cells) and adaptive immune cells (T cells and B cells) play an important role in tumor progression (Thakkar et al., 2020). High CD10 + /low CD20 + immune cell infiltration ratio is an important prognostic factor for lung squamous cell carcinoma (Kadota et al., 2015). Immunomodulatory therapy of tumor-specific neutrophil and B lymphocyte responses may be suitable in the treatment of lung squamous cell carcinoma (Kadota et al., 2015) (65). Correale P' team investigated the prognostic value of tumor-infiltrating CD8 + T cells expressing the chemokine receptor 7 [T(ccr7)], demonstrating that patients with colorectal cancer with high T (ccr7) and T (reg) invasion have a better prognosis (Correale et al., 2012). In the present study, the expression of prognostic CRGs in HCC was significantly and positively correlated with the infiltration of immune cells. Therefore, the target of prognostic CRGs with the aim of interfering with the function of immune cell infiltration in HCC might provide a new solution for immunotherapy against HCC.
At present, immunotherapy targeting immune checkpoints, especially programmed cell death protein 1/programmed cell death ligand 1 (PD-1/PD-L1) and CTLA-4 blockers have become feasible in the treatment of many malignant tumors (Cai et al., 2019). PD-1 is a member of the CD28 family. PD-1 and its ligands are widely expressed in T cells and play a broader immunoregulatory role in T cell activation and tolerance (Cha et al., 2019). PD-L1 is a transmembrane protein considered as a co-suppressor of the immune response, and it binds to PD-1 to reduce the proliferation of PD-1-positive cells, inhibits their cytokine secretion and induce apoptosis. PD-L1 also plays an important role in various malignancies, attenuating the host immune response to tumor cells (Han et al., 2020). The PD-1/PD-L1 axis is responsible for cancer immune escape and has a huge impact on cancer therapy, and its inhibition is an effective treatment for many cancers (Dermani et al., 2019). Previous studies found that copper regulates a key signaling pathway that mediates PD-L1-driven cancer immune evasion (Voli et al., 2020). Cytotoxic T lymphocyte-associated antigen 4 (CTLA-4) is a membrane glycoprotein expressed by activated effector T cells and is involved in the inhibition of T cell proliferation, cell cycle progression and cytokine production (Zhao et al., 2018).
Antibodies targeting CTLA-4 or in combination with other therapies significantly enhance the antitumor effects and improve the prognosis in malignant diseases (Zhao et al., 2018). This study found that the expression of prognostic CRGs in the TCGA-LIHC cohort was positively correlated with the expression of PD-1 (PDCD1), PD-L1 (CD-274) and CTLA-4, suggesting that the prognosis of HCC patients with high expression of CRGs should be improved by immunotherapy targeting immune checkpoints.
The research on tumor immune checkpoint inhibitors in tumor immunotherapy is the most developed. TMB and MSI are considered potential predictive biomarkers involved in response to ICIs (Dudley et al., 2016;Ritterhouse, 2019). High TMB is associated with the response to ICI in multiple tumor types (Weir et al., 2004). Our results showed that CDKN2A increased the TMB score in HCC and PDHA1 increased the MSI score. Furthermore, prognostic CRGs were positively or negatively associated with multiple chemotherapeutic agents. These results provide new potential therapeutic targets in the treatment of HCC.
N6-methyladenosine (m6A) is the most prevalent internal mRNA modification in mammalian cells. RNA methylation is a pervasive post-transcriptional modification that plays a key role in regulating various biological processes, and its dysregulation is closely related to the occurrence of human malignancies (Anita et al., 2020;Chen and Wong, 2020;Han et al., 2021;Liu et al., 2021;Zhang et al., 2021) through various mechanisms, providing more possibilities for the early diagnosis and treatment of cancer (Sun et al., 2019). METTL3 is involved in pancreatic carcinogenesis and is a potential prognostic marker or therapeutic target (Xia et al., 2019). HBXIP promotes the progression of gastric cancer via METTL3mediated MYC mRNA m6A modification . Currently, the link between CRGs and m6A-related genes in HCC has not been investigated, which was one of the aims of the present study. Our result showed that m6A-related genes were significantly upregulated in the C1 cluster with higher expression of CRGs, and prognostic CRGs were significantly correlated with m6A-related genes, suggesting that CRGs might affect the progression of HCC through m6A modification. However, further studies should be performed to confirm this result.
Stromal components, including CAFs and TAMs, play important roles in cancer initiation and progression (Hanahan and Coussens, 2012). Our results showed that prognostic CRGs upregulated the expression of macrophages, fibroblasts, endothelial cells, and hepatic progenitor cells in HCC. Furthermore, prognostic CRGs were positively correlated with many markers of CAFs and macrophages. Previous studies showed that CAFs promote the progression of HCC (Ju et al., 2009;Affo et al., 2021). Clinical studies and experimental mouse models also strongly suggest that TAMs promote tumor progression (Qian and Pollard, 2010;Noy and Pollard, 2014). The present study found that prognostic CRGs upregulated the expression of macrophages, fibroblasts, endothelial cells, and hepatic progenitor cells in HCC. Moreover, prognostic CRGs were significantly and positively correlated with many markers of stromal cells, especially CAFs and TAMs. Therefore, prognostic CRGs might affect the progression of HCC patients by altering the expression of CAFs and TAMs in the TME.
Eight lncRNA-miRNA-mRNA networks were also constructed, lncRNA DANCR reduces the expression of miR-125b-5p and regulates the proliferation and apoptosis of hepatoma cells (Yang Frontiers (Bo et al., 2021). LncRNA DLEU2 aggravates the progression of HCC . Furthermore, the GSEC/ miR-101-3p/SNX16/PAPOLG axis of the ceRNA network axis is an important factor associated with HCC prognosis and immune infiltration (Hu et al., 2022). Our study revealed that these mRNA-miRNA-lncRNA networks were associated with the prognosis of HCC patients. All these pieces of evidence suggest that these regulatory axes might play an important role in the progression of HCC. However, further studies should be performed to confirm this result.

Conclusion
CRGs affected the infiltration of immune cells such as macrophages and CD8 + T cells, promoted the expression of immune checkpoints such as CD274 and HAVCR2, and promoted the expression of various m6A-related genes in HCC. CRGs also participated in multiple cancer-related signaling pathways in HCC. A total of 19 CRGs that were highly expressed in HCC were analyzed, and five CRGs (CDKN2A, DLAT, DLST, GLS, and PDHA1) that were associated with the prognosis of HCC patients were identified. A prognostic model of CRGs was established, which predicted more accurately OS and PFS in HCC patients. Mutations in prognostic CRGs led to poorer prognosis in patients with HCC. Prognostic CRGs were positively correlated with B cells, T cells, macrophages, and other immune cells, and were positively correlated with various chemotherapeutic drugs. They affected the TME in HCC by affecting CAFs and TAMs. Finally, eight lncRNA-miRNA-mRNA regulatory axes that affected the progression of HCC were predicted. Hence, our study laid the foundation for an in-depth understanding of the role of CRGs in HCC.

Author contributions
Conceptualization, HZ, JW and FY; methodology, QY and YZ; software, GZ and SY; formal analysis, JW and YY; investigation, HQ, WS and QW; data curation, YC and PZ; writing-original draft preparation, HQ, WS, GZ; writing-review and editing, FY, JY; supervision, FY, JW, JY; funding acquisition, FY, JY. All authors have read and agreed to the published version of the manuscript.